//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 1998-2002 by Jeongnim Kim, Lars Jonsson and Seungwon Lee
//   All rights reserved.
//////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-

#ifndef OHMMS_UBLAS_ARPACKSYMENGINE_H
#define OHMMS_UBLAS_ARPACKSYMENGINE_H

#include <complex>
#include <iostream>
using namespace std;

#include "EigSolver/EigenSpaceDimension.h"
#include "EigSolver/Arpackf.h"
#include "OhmmsNumeric/OhmmsBlasf.h"
#include <boost/numeric/ublas/config.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/traits.hpp>

template<class T>
class BoostArpackSymSolver {

public:

  typedef T value_type;
  typedef boost::numeric::ublas::upper                              upper;
  typedef typename boost::numeric::ublas::type_traits<T>::real_type real_type;
  typedef boost::numeric::ublas::matrix<value_type>                 EigenStates_t;
  typedef boost::numeric::ublas::unbounded_array<value_type>        WorkArray_t;
  typedef boost::numeric::ublas::unbounded_array<real_type>         RealArray_t;
  typedef boost::numeric::ublas::unbounded_array<int>               IndexArray_t;
  typedef typename WorkArray_t::size_type                           size_type;

  EigenStates_t Psi;

protected:

  EigenSpaceDimension& States;
  int N;
  int M;
  int ncv; // number of Lanczos vectors
  int info;
  int lworkl;
  int ido;
  int ierr;
  int nconv;
  int maxiter;
  int ishfts;
  int mode;
  int rvec;
  string which;
  int iparam[11], ipntr[14];

  IndexArray_t select;
  WorkArray_t workd; // array of 3N+1
  WorkArray_t workl; // array of lworkl
  WorkArray_t resid; // residue
  WorkArray_t eigVal;// eigen values
  WorkArray_t eigVec;// eigen vectors 
  RealArray_t rwork; // only for a complex problem

  real_type tol;  // tolerance
  value_type sigma; // reference value for shift-invert

public:


  BoostArpackSymSolver(EigenSpaceDimension& dim, bool nonortho=false):
    States(dim) {   
    tol = 1.0e-13;
    info = 0;
    ido = 0;
    maxiter = 1000;
    ncv = -1;
  }

  ~BoostArpackSymSolver() { }

  template<class HamiltonianMatrix_t>
  bool getEigvectors(HamiltonianMatrix_t& H) {

    initializeWorkSpace();
    do{
      solve(eigVec.begin());
      H.apply(workd.begin()+ipntr[0]-1, 
	      workd.begin()+ipntr[1]-1);
    } while(ido == -1 || ido == 1);
    rvec = 1;
    finalize(eigVec.begin());

    if(info < 0) {
      cout << endl << "Error with __aupd, info = " <<  info << endl;
      cout << " Check documentation in __aupd " << endl << endl;
      return false;
    } else
      return true;
  }

  template<class HamiltonianMatrix_t>
  bool getEigvalues(HamiltonianMatrix_t& H) {

    initializeWorkSpace();
    do{
      solve(eigVec.begin());
      H.apply(workd.begin()+ipntr[0]-1, 
	      workd.begin()+ipntr[1]-1);
    } while(ido == -1 || ido == 1);

    rvec = 0;

    finalize(eigVec.begin());

    if(info < 0) {
      cout << endl << "Error with __aupd, info = " <<  info << endl;
      cout << " Check documentation in __aupd " << endl << endl;
      return false;
    } else
      return true;
  }

  /////////////////////////////
  // currently not all arpack functions are implemented
  /////////////////////////////
  template<class HamiltonianMatrix_t>
  bool getEigvalues(HamiltonianMatrix_t& H, HamiltonianMatrix_t& S) {
    return false;
  }

  template<class HamiltonianMatrix_t>
  bool getEigvectors(HamiltonianMatrix_t& H, HamiltonianMatrix_t& S) {
    return false;
  }

  template<class HamiltonianMatrix_t, class IterType>
  bool getEigvalues(HamiltonianMatrix_t& H, IterType phi) {
    return false;
  }

  template<class HamiltonianMatrix_t, class IterType>
  bool getEigvectors(HamiltonianMatrix_t& H, IterType phi) {
    return false;
  }

  template<class HamiltonianMatrix_t, class IterType>
  bool getEigvalues(HamiltonianMatrix_t& H, HamiltonianMatrix_t& S, 
		    IterType phi) {
    return false;
  }

  template<class HamiltonianMatrix_t, class IterType>
  bool getEigvectors(HamiltonianMatrix_t& H, HamiltonianMatrix_t& S,
		     IterType phi) {
    return false;
  }


  template<class HamiltonianMatrix_t>
  bool checkEigvectors(HamiltonianMatrix_t& H) {
    if ( ierr != 0) {
      cout << endl << "Error with __eupd, info = " << ierr << endl;
      cout << "Check the documentation of __eupd. " << endl;
      return false;
    } else {
      
      nconv =  iparam[4];
      if(rvec == 1) {
	WorkArray_t ax(N);
	WorkArray_t res(M);
	
	for(int j=0; j<nconv; j++) { // residual calculations
	  H.apply(eigVec.begin()+j*N, ax.begin());
	  value_type da  = -eigVal[j];
	  int incx = 1;
	  int incy = 1;
	  OhmmsBlas::axpy(N, da, eigVec.begin()+j*N, incx, ax.begin(), incy);  
	  res[j] = OhmmsBlas::norm2(N, ax.begin(), incx);
	  res[j] /= fabs(eigVal[j]);
	}
	cout << "Ritz values and relative residuals " << endl;
	for(int j=0; j<nconv; j++)
	  cout << j << "  "  << eigVal[j] << "  " << res[j] << endl;
      } else {
	cout << "Ritz values only " << endl;
	for(int j=0; j<nconv; j++)
	  cout << j << "  "  << eigVal[j] <<  endl;
      }

      cout << endl;
      if(info == 1) {
	cout << endl << " Maximum number of iterations reached." << endl;
      } else if ( info == 3) {
	cout << endl << "No shifts could be applied during implicit "
	     <<   "  Arnoldi update, try increasing NCV." << endl;
      }
      cout << endl;
      cout << " _SDRV1 " << endl;
      cout << " ====== " << endl << endl;
      cout << " Size of the matrix is " << N << endl;
      cout << " The number of Ritz values requested is " << M << endl;
      cout << " The number of Arnoldi vectors generated  (NCV) is " << ncv << endl;
      cout << " What portion of the spectrum: " << which << endl;
      cout << " The number of converged Ritz values is " << nconv << endl;
      
      cout << " The number of Implicit Arnoldi update iterations taken is ";
      cout <<  iparam[2] << endl;
      cout << " The number of OP*x is " << iparam[8] << endl;
      cout << " The convergence criterion is " << tol << endl;
      return true;
    }
  }

  inline void setMaxIter(int n) { maxiter = n;}
  inline void setMaxTrialStates(int n) { ncv = n;}
  inline void setTolerance(real_type x) { tol = x;}
  inline real_type eigval(size_type i) const { return eigVal[i];}

private:
  void initializeWorkSpace() {
    // getting dimension from EigenSpace
    N = States.N;
    M = States.M;

    if(ncv < 0) ncv = 4*M;
    if(ncv > N) ncv = N;
    //TODO: COMPLEX and DOUBLE  (maybe not)
    lworkl = ncv*(3*ncv+5);
    iparam[0] = 1; iparam[6] = 1;
    which = "SM";
    iparam[2] = maxiter;
    if(workd.size() < 3*N) {
      workd.resize(3*N);
    }                                    
    if(resid.size() < N) {
      resid.resize(N);
    }
    if(select.size() < ncv) {
      select.resize(ncv);
    }
    if(eigVal.size() < ncv*2) {
      eigVal.resize(ncv*2);
    }
    if(eigVec.size() < ncv*N) {
      eigVec.resize(ncv*N);
    }
    if(workl.size() < lworkl) {
      workl.resize(lworkl);
    }
    
  }

  void solve(double* eigv){
    int n = N; int ldv = N; int nev = M;
    char bmat = 'I';
    dsaupd(&ido, &bmat, &n, which.c_str(), &nev, &tol,
           resid.begin(), &ncv, eigv, &ldv, iparam, ipntr,
           workd.begin(), workl.begin(), &lworkl, &info);
  }

  void solve(complex<double>* eigv){

    if(rwork.empty()) {
      rwork = RArray_t(ncv); // only for a complex problem
    } 
    int n = N; int ldv = N; int nev = M;
    char bmat = 'I';
    znaupd(&ido, &bmat, &n, which.c_str(), &nev, &tol,
           resid.begin(), &ncv, eigv, &ldv, iparam, ipntr,
           workd.begin(), workl.begin(), &lworkl, rwork.begin(), &info);
  }

  void solve(float* eigv){
    int n = N; int ldv = N; int nev = M;
    char bmat = 'I';
    ssaupd(&ido, &bmat, &n, which.c_str(), &nev, &tol,
           resid.begin(), &ncv, eigv, &ldv, iparam, ipntr,
           workd.begin(), workl.begin(), &lworkl, &info);
  }

  void finalize(double *eigv) {
    char bmat = 'I';
    char howmany ='A';
    int n = N; int ldv = N; int nev = M;
    dseupd( &rvec, &howmany, select.begin(),
            eigVal.begin(), eigv, &ldv, &sigma,
            &bmat, &n, which.c_str(), &nev, &tol, resid.begin(), &ncv,
            eigv, &ldv, iparam, ipntr,
            workd.begin(), workl.begin(), &lworkl, &ierr );

    Psi.resize(N,M);
    OhmmsBlas::copy(N*M,eigv,Psi.data().begin());
    //for(int i=0; i<States.M; i++) ene[i] = eigVal[i];
  }

  void finalize(complex<double> *eigv) {
    char bmat = 'I';
    char howmany ='A';
    int n = N; int ldv = N; int nev = M;
    WorkArray_t workev(2*ncv);    
    zneupd( &rvec, &howmany, select.begin(),
            eigVal.begin(), eigv, &ldv, &sigma, 
            workev.begin(),  &bmat, &n, which.c_str(), &nev, &tol, 
            resid.begin(), &ncv, eigv, &ldv, iparam, ipntr, workd.begin(), 
            workl.begin(), &lworkl, rwork.begin(), &ierr );
    //for(int i=0; i<States.M; i++) ene[i] = eigVal[i].real();
  }

  void finalize(float *eigv) {
    char bmat = 'I';
    char howmany ='A';
    int n = N; int ldv = N; int nev = M;
    sseupd( &rvec, &howmany, select.begin(),
            eigVal.begin(), eigv, &ldv, &sigma,
            &bmat, &n, which.c_str(), &nev, &tol, resid.begin(), &ncv,
            eigv, &ldv, iparam, ipntr,
            workd.begin(), workl.begin(), &lworkl, &ierr );
    //for(int i=0; i<States.M; i++) ene[i] = eigVal[i];
  }
};


#endif

/***************************************************************************
 * $RCSfile: BoostArpackSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: BoostArpackSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

